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Abstract 

In this paper, we propose a variational formulation to study the singular evolu- 
tion equations that govern the dynamics of surface modulations on crystals below 
the roughening temperature. The basic idea of the formulation is to expand the 
surface shape in terms of a complete set of basis functions and to use a varia- 
tional principle equivalent to the continuum evolution equations to obtain coupled 
nonlinear ordinary differential equations for the expansion coefficients. Unlike sev- 
eral earlier approaches that rely on ad hoc regularization procedures to handle the 
singularities in the evolution equations, the only inputs required in the present ap- 
proach are the orientation dependent surface energies and the diffusion constants. 
The method is applied to study the morphological equilibration of patterned unidi- 
rectional and bidirectional sinusoidal modulations through surface diffusion. In the 
case of bidirectional modulations, particular attention is given to the analysis of the 
profile decay as a function of ratio of the modulating wavelengths in the coordinate 
directions. A key question that we resolve is whether the one-dimensional decay 
behavior is recovered as one of the modulating wavelengths of the two-dimensional 
profiles diverges, or whether one-dimensional decay has qualitatively distinct fea- 
tures that cannot be described as a limiting case of the two-dimensional behavior. 
In contrast to some earlier suggestions, our analytical and numerical studies clearly 
show that the former situation is true; we find that the one-dimensional profiles, 
like the highly elongated two-dimensional profiles, decay with formation of facets. 
While our results for the morphological equilibration of symmetric one-dimensional 
profiles are in agreement with the free-boundary formulation of Spohn, the present 
approach can also be used to study the evolution of asymmetric profile shapes 
where the free-boundary approach is difficult to apply. The variational method 
is also used to analyze the decay of unidirectional modulations in the presence of 
steps that arise in most experimental studies due to a small misorientation from 
the singular surface. 

Keywords: Models of nonlinear phenomena. Surface diffusion. Stepped single 
crystal surfaces 
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1 Introduction 



In the recent years, nanoscale material structures have been fabricated on soHd surfaces using 
various patterning techniques which, in some cases, exploit the phenomenon of stress-driven 
self-assembly. Potential application of these structures as functional elements in nano- and 
micro-electronics has led to considerable interest in understanding their stability and the 
processes that play a role in determining their morphological evolution. Due to their small 
size, the time scales involved in the processing and fabrication of these structures are small. 
Their small size also implies that the surface energy is a significant part of the overall energy 
of individual structures and, hence, is important in determining their morphological evolution. 
For example, transport of mass via surface diffusion can lead to smaller surface areas, thereby 
reducing surface energy for orientation independent surface energy density. Surface diffusion 
can also result in transformation of surface orientations with large surface energies to certain 
preferred low-energy orientations. Therefore, a fundamental understanding of the kinetics of 
surface transport, including the role of orientation-dependent surface energy, is important in 
the fabrication of nanostructures. 

In the case of patterned features on the surfaces of amorphous solids or on crystalline 
surfaces above the thermodynamic roughening temperature, the surface energy is independent 
of the surface orientation. In either situation, surface diffusion is driven solely by the ten- 
dency to reduce the surface area and is described within the framework established through 
the classic work of Herring |Q and Mullins |^ . According to their theory, a small sinusoidal 
perturbation on a flat surface retains its shape while the amplitude decays at a rate that is 
determined by the wavelength of the perturbation. In the case of a crystal below its ther- 
modynamic roughening temperature, on the other hand, the situation is more complicated. 
Because a surface with an orientation that is close to a high symmetry orientation is made 
up of steps, the dependence of surface energy density on orientation acquires a cusp at such 
a high symmetry orientation. 

Experiments have shown that sinusoidal perturbations about these high symmetry sur- 
face orientations do not retain the sinusoidal shape. Instead, the peaks in the profile become 
flattened, with the extent of the flat tops increasing as the amplitude of the profile decays 
||3-10|. The strong nonlinearity associated with the cusp in the surface energy causes difficul- 
ties in extending the theory of Mullins to faceted crystals. In this case, the decay of periodic 
surface profiles has been theoretically studied using various continuum models p[j5,1116|, by 
explicitly accounting for the dynamics of individual steps []l^|l7|jT^ and by using Monte-Carlo 
simulations |19-2^. While useful insights have been gained, each method has its shortcomings 
and the problem of decay remains open |2|-28|. In order to motivate the present work, which 
is based on a variational framework for handling the nonlinear and singular aspects of surface 
evolution of crystals below the roughening temperature, we briefiy summarize here the key 
results obtained in these studies. 

The cusp in the surface energy leads to partial differential equations for surface evolution 
that involve Dirac-delta functions or other nonanalytic functions of the surface slope. In the 
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case of both unidirectional and bidirectional sinusoidal perturbations, these equations have 
been solved by smoothing out the singularities |^^JTl|. Even though calculations performed 
using this procedure show flattening of profile shapes due to facet formation, in agreement 
experimental observations, the approach is ad hoc and it is not clear that different smoothing 
schemes would lead to identical results. For unidirectional modulations, an alternate approach 
was taken by Villain, Zangwill and coworkers 1 12 -141, who argued that the surface chemical 
potential for surface diffusion should not involve the singular term that arises from the cusp. 
Retaining only the nonsingular contributions arising from step interactions, they showed that 
the sinusoidal profile does not flatten as it decays but, to the contrary, the top of the profile 
becomes sharper with time. According to their model, the curvature of the profile diverges 
at its extremum points. In distinct contrast to the above approaches, Spohn developed an 
ingenious method to deal with the singular nature of the evolution equations for unidirectional 
modulations by transforming them to a free-boundary problem [IT^Jl^ . In his approach, if 
the location of the facet is known in advance (which is true for highly symmetric profiles), an 
ordinary differential equation for the length of the facet can be coupled to nonsingular partial 
differential equations for the evolution of stepped regions to obtain a coupled problem that is 
well-posed. Numerical solution of these equations in the case of unidirectional profiles shows 
pronounced flattening of the tops and bottoms of the periodic profiles. 

In order to resolve the inconsistent predictions of the various continuum theories, more 
atomistic based approaches, particularly Monte Carlo (MC) simulations |19-2^ and models 
that incorporate the dynamics of straight steps |14,1^, have been pursued by several investi- 
gators. While some of the MC simulations show greater tendency to formation of facets than 
do others, none of the simulations indicate a sharpening of the profile as predicted in | ]l^ 
and Since kinetics of surface diffusion within the MC simulations is slow, most of these 
simulations have been confined to small sample sizes which are sensitive to step fluctuations. 
Furthermore, most of the MC simulations do not account for the long-range interactions of 
surface defects and, therefore, cannot capture the inverse-square elastic or electrostatic in- 
teractions between steps. This problem does not occur in step dynamics models where step 
interactions can be included explicitly. However, this approach also has limitations, since 
the decay kinetics and profile shapes depend strongly on the details of the short-range step 
interactions which are poorly understood. For example, if an attractive interaction is present 
between unlike steps at the top of the sinusoidal profile, the profile tends to be flatter and 
decays more rapidly @|. Also, in the case of bidirectional profiles, step dynamics models 
become very complicated and tractable solutions have been achieved only for cylindrically 
symmetric surface perturbations that involve circular steps |17|; only approximate analytic 
results for the case of surface evolution due to evaporation kinetics have been obtained for 
bidirectional sinusoidal perturbations [l^. The more complicated and experimentally rele- 
vant case of surface diffusion limited decay of sinusoidal profiles has not been treated in a 
quantitative way using discrete step-dynamics models. 

In this paper, we propose a variational approach to solve the nonlinear evolution equa- 
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tions for surface evolution and address the problem of decay of both one- and two-dimensional 
sinusoidal profiles on crystals below the roughening temperature. In the past, variational ap- 
proaches have been used to study surface evolution involving complex geometries and shapes, 
but attention was limited to functional forms of surface energy that are well-behaved. For 
example, a variational model for evolution of a surface with an isotropic surface energy was 
introduced by Needleman and Rice to study the problem of grain boundary transport 
and cavitation. In recent years, Suo has applied the variational formulation to con- 
sider surface evolution during motion of defects such as cracks, voids and inclusions. The 
formulation of the variational approach adopted here is similar, and attention is focused on 
addressing the difficulties posed by nonanalytic forms of the surface energy. The basic idea of 
the formulation is to expand the surface shape in terms of a complete set of basis functions 
and to use a variational principle derived from the continuum evolution equations to obtain 
coupled nonlinear ordinary differential equations for the expansion coefficients. The advan- 
tage of this method is that even though the evolution equations are singular, the ordinary 
differential equations that determine the expansion coefficients are well-behaved and can be 
integrated by standard numerical techniques. The only inputs required to solve them are the 
orientation-dependent surface energies and the surface diffusion constants; there are no other 
ad hoc parameters in the model. The conclusions drawn through application of this approach 
are briefly summarized here and the details are described in the sections to follow. 

In general, the main strength of the variational approach is that it provides a basis for 
extracting approximate solutions of a set of equations representing a physical phenomenon 
without the need for a priori choices on which features of those equations to retain or discard. 
The approximation made in this approach is usually the introduction of a lower limit on 
resolution of features of the solution, so the approximate solution found is the "best" available 
within this restricted range of possible solutions. The lower limit on resolution can be reduced 
indefinitely, in principle, and convergence is expected. However, convergence can be proved 
in only the simplest circumstances, and we are normally content with establishing essential 
features of behavior, with the physics of the system providing guidance on the level of detail 
needed. As will be demonstrated below, extensive numerical testing with ever finer resolution 
is a good indicator of the level of detail required in any particular problem. 

In the case of decay of unidirectional sinusoidal profiles, when the step creation energy 
is finite, our results agree with the predictions of Spohn [15,^]. Here, we find that the decay 
rate of the amplitude increases with the magnitude of step formation energy. On setting the 
step creation energy to zero (or if the singular term in the chemical potential is ignored), 
we find perfect agreement with the results of Villain and Zangwill |1^,14|. This implies 
that the term arising from the cusp in the surface energy represents a singular perturbation; 
the profile decay in the presence of the this term, however small, is qualitatively different 
from the decay when this term is absent. While our model is in agreement with the free- 
boundary approach for symmetric profiles, it has several advantages from both computational 
and conceptual viewpoints. For example, the free boundary approach requires the knowledge 
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of the location of facets at the beginning of a calculation, but this is not known for the case 
of asymmetric profiles. Within the variational formulation, this information is not needed a 
priori and the approach can therefore be applied to study the evolution of arbitrary profile 
shapes. The numerical procedure for solving the free-boundary problem involves both partial 
and ordinary differential equations, and these solutions must be matched through certain 
continuity conditions at the edges of the facet. The present variational formalism involves 
only solution of coupled ordinary differential equations and is therefore more straightforward 
in its implementation. Perhaps the most significant advantage of the variational approach is in 
handling the evolution of two-dimensional profiles, where the formulation of the free-boundary 
problem becomes mathematically very involved. We are not aware of any calculations of two- 
dimensional decay of sinusoidal profiles, other than by regularizing the nonanalytic features 
in the surface energy. Once again, the present approach does not require ad hoc parameters 
and is not fundamentally different from the one-dimensional case. 

The ability to handle two-dimensional profiles allows us to address an important and 
outstanding issue concerning the limiting behavior of highly elongated two dimensional pro- 
files, that is, profiles with a wavelength of modulation in one of the coordinate directions that 
is much larger than the wavelength in the other direction. The question we would like to 
answer is whether the one-dimensional decay behavior is recovered as one of the modulat- 
ing wavelengths of the two-dimensional profiles diverges, or whether one-dimensional decay 
has qualitatively distinct features that cannot be described as a limiting case of the two- 
dimensional behavior. It has been suggested by Rettori and Villian |12| that the later case is 
true, who point out the facets form in the case of two dimensional profiles because of the line 
tension contribution of steps in the surface chemical potential. Since there is no contribution 
of the step formation energy in their chemical potential for the ID profiles, they conclude 
that there is no facet formation in the ID case. Within the present approach, where the 
line tension contribution of the steps is explicitly accounted for, we are able to analytically 
establish that the limiting behavior of the 2D profiles approach the ID case. In particular, 
we show that the evolution equations of the variational parameters that express the shape of 
the 2D profiles approach the corresponding equations for the ID profiles as the aspect ratio 
of the 2D profiles (the ratio of wavelengths in the coordinate directions) approaches infinity. 
Our analytical results are confirmed by numerical calculations which show that the 2D profile 
decay is already close to the ID behavior when the aspect ratio is between 100 and 1000. 

In most experimental situations, the study of decay of the ID sinusoidal perturbations 
is infiuenced by the presence of widely spaced steps that run in a direction perpendicular to 
the direction of the corrugation. These steps arise as a result of a very small miscut from the 
singular orientation and are difficult to eliminate]^. In the present work we will show that, 
as the miscut angle tends to zero, the decay profiles approach the shape of the sinusoidal ID 
profile on the singular surface. In contrast to earlier expectations ||l^, our studies indicate 
that the tendency to form facets persists even as the miscut slope goes to zero. 



^STM images that show the presence of these steps on patterned Si and Au surfaces are given in 
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This paper is organized in the following away. Sec. 2 is devoted to the development of 
the variational formulation for both one and two-dimensional profiles. Here, we will derive 
evolution equations in the case where mass transport is controlled by terrace-diffusion. Ana- 
lytical calculations establishing the limiting one-dimensional behavior of highly elongated two 
dimensional profiles will also be presented in this section. In Sec. 3, we will present numerical 
results for the decay of one-dimensional and two-dimensional profiles and profiles on surfaces 
that are miscut with respect to a singular surface. Particular emphasis is given to the anal- 
ysis of the decay rate of two-dimensional profiles as a function of their aspect ratios and the 
behavior of sinusoidal profiles on miscut surfaces as the angle of miscut vanishes. In the case 
of one-dimensional profiles, we will also consider the decay of asymmetric profile shapes and 
present numerical results on the convergence properties of the variational formulation. Some 
the predictions of our calculations are compared with recent experiments [31|. A summary of 
key results is given in Sec. 5. 



2 Variational Formulation 

We start by considering the free energy F of the surface S, written in terms of the orientation- 
dependent surface energy 7 as 

F = / 7 (V5/1) dS, (1) 
Js 

where h denotes the height of the surface measured relative to a high symmetry reference 
plane and Vg' is the interior gradient operator on the surface. Using Vn to denote the normal 
velocity of the surface relative to the material points instantaneously on that surface, the rate 
of change of the total free energy can be expressed in terms of the surface chemical potential 
H as 

F = [ f,VndS. (2) 
Js 

The above equation provides a definition for the surface chemical potential, which is a surface 
field representing the driving force for change of shape of the crystal by diffusive mass transport 
over its surface. The chemical potential can be viewed as a continuum representation of the 
sensitivity of the free energy of the system to alterations in surface shape. The mass flux on 
the surface, here denoted by the surface vector field j, is assumed to be related to the surface 
chemical potential through a kinetic relationship of the form 

j = -CVs/i, (3) 

where C is a positive mobility parameter that, among other things, can depend on the local 
surface slope, temperature and the concentration of the diffusing species.^ Conservation of 

^Here, for convenience, we assume that surface diffusion is isotropic, so that the kinetic relationship is 
written using only one mobility parameter. The variational formulation that will be developed in this section 
can be generalized to include diffusion anisotropy in a straightforward manner. 
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mass at each point on the surface requires that the net rate of accumulation of mass at that 
point on the surface, determined by the normal speed is the negative of the divergence of 
the surface mass flux at that point, or 

Vn = -^S-j- (4) 



When the surface energy is sufficiently smooth, eqns can be solved by standard numer- 
ical techniques to obtain the evolution of the surface shape. Alternatively, the weak form of 
the equations can be solved by employing a finite-element formulation or other variational 
approaches. However, if the dependence of a surface energy on surface shape is not ana- 
lytic, a direct solution of the evolution equations becomes difficult due to the singularities 
in the chemical potential. In such cases, variational approaches provide a powerful tool for 
determining the evolution of the surface. 

The development of the variational formulation involves two key ingredients. First, 
suppose that q is any possible mass flux field, sharing properties of periodicity and smoothness 
with j, but is otherwise arbitrary. If we take the inner product of each side of (^) with q and 
integrate each side over S, after the application of a vector identity, the result can be written 
in the form 

/ Vs ■ {m)dS - [ fiVs- qdS + [ ^j-qdS = 0. (5) 

Js Js Js 

If the fields and the surface shape are spatially periodic and if S includes exactly one period 
of the system, then the first term in the above equation vanishes. This follows from the fact 
that fi would have the same value on opposite sides of a segment of S due to periodicity, and 
the mass flux out of S would exactly cancel mass flux out on the opposite side for the same 
reason. This result provides the basis for an extremum principle which can be expressed as a 
maximum principle or a minimum principle; the two forms are equivalent and we adopt the 
latter form here. 

Next, define a functional over the range of admissible trial fields j* as 



dS. (6) 



Then (|5|) implies that the functional is stationary under arbitrary variations of the surface flux 
j* when the flux field is the actual flux field, that is, when j* = j. This follows immediately 
from (^) if q is interpreted as 6j* in the usual variational notation. The arbitrariness of (5j* 
and the fundamental theorem of the calculus of variations require that j must satisfy (^) 
pointwise on 5*. The fact that $[j*] is not only stationary but also an absolute minimum 
under variations about j can be seen by considering $[j* -|- — to second order in 5j*. 
Since the mobility parameter C is positive, it follows from @ that 

m* + m - m*] = ■ srds > o, (?) 

thus ensuring that j* = j is indeed a minimum of 
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The main advantage of the variational formulation is that instead of solving the con- 
ventional evolution equations, we can expand the mass flux in terms of a complete set of basis 
functions and obtain the time dependence of the expansion coefficients by minimizing the 
functional There is no need to deal explicitly with the singular chemical potential while 
carrying out this minimization process. This can be seen by considering the first term in (^) 
which can be seen to be equal to F if we replace the normal velocity f„ in favor of — V5 • j* 
in (|^; F can be directly evaluated in terms of j* without using the chemical potential. The 
variational functional can then be conveniently expressed as 

$[j1 = i^[j1+ l^i*-i*ds. (8) 

Below, we provide details of the implementation of the variational formulation to study the 
decay of periodic one and two dimensional profiles; the one-dimensional case is considered 
first. 



2.1 One-dimensional profiles 

The orientation-dependent surface energy that applies to one-dimensional modulations can 
be written as 

l{hx) = Pi\K\ + ^\h,f, (9) 

where the subscript x denotes the derivative with respect to the spatial coordinate and the 
terms proportional to Pi and P3 represent step formation and step interaction energies, re- 
spectively. This form of the surface energy expression presumes that \hx\ ^ 1. If we focus 
attention on periodic profiles with wavelength A, then free energy in one period and the 
position dependent surface chemical potential can be written as 

F= (^-i{K)dx, (10) 







and ^ ^ 

= ~~r ( TT ) = ~(^^^ [^J ~ 2/33 \hx\ hxx , (11) 



dx \dhx, 

respectively, where 6[ ] denotes the Dirac delta function. The singular term in the chemical 
potential can be ignored if the curvature vanishes at the points where the slope goes to zero. 
However, since the step interaction term is cubic in the surface slope, it can be shown that 
the curvature diverges like {x — xo)~^^^ around an extremum located at x = xq |12.14]. This 
implies that the first term plays a crucial role in determining the decay shape. 

If the surface shape is symmetric about the origin, it can be expanded in terms of a 
Fourier cosine series as 

N 

h{x, t) = ^ a„(t) cos{nkx), < x < A, (12) 
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where k = 2-k/\ and the number of Fourier modes is finite but otherwise unrestricted.^ 
If attention is focused on small perturbations that satisfy the criterion \hx\ <C 1, the normal 
velocity of the surface can be written in terms of surface shape as Vn = ht, where subscript t 
denotes differentiation with respect to time. We can now invoke conservation of mass (Q) to 
express the surface mass flux as 

./ N - / .sm(nkx) 

j{x,t) = -Y,an{t) \ ' . (13) 

n=l 

The variational formulation provides first-order coupled ordinary differential equations for the 
Fourier coefficients an{t); these equations will be numerically integrated to obtain the evolving 
surface shape for t > using information on the surface shape at t = 0. 

In order to write the functional <I> (refer to (^) in terms of the variational parameters, 
which in the present case are the dnS, we first have to specify the dependence of the surface 
mobility parameter on the surface slope. This dependence can be obtained by considering 
the mechanisms involved in determining mass transport on stepped surfaces, namely, the 
attachment / detachment of atoms from step-edges and the diffusion of adatoms on the terraces 
separating the steps. In this paper we will restrict attention to the case where the mass 
transport is limited by terrace diffusion; the mobility parameter in this situation does not 
depend the surface slope.Q 

Using the expression for mass flux in (H), we find the variational functional for the 
terrace-diffusion limited case to be 

N N -2 

$[di,. . . ,(iAr] = - ^ Hn[ai, . . . ,Oiv]d„(t) + —3— ^ (14) 



n=l n=l 



where 



Hn[ai, . . . ,aN] = nkPi / Sgn [h^] sm{nkx)dx + nkfd-s / Sgn[hx]\hx\^ sm{nkx)dx, (15) 

Jo Jo 

in which Sgn[p] = d\p\/dp. Minimizing the functional by setting d^/ddn = for n = 1, . . . , A^, 
we obtain the evolution equations 

dn{t) = ^ ^ Hn[ai, . . . , ajv] (16) 
vr 

for the Fourier coefficients in ([l^). If the initial shape is known, approximate solutions to 
these coupled first-order ordinary differential equations can be found using standard numerical 
integration techniques. 

^Numerical studies on the sensitivity of the results to the number of Fourier coefficients and the convergence 
properties of the present formulation will be presented in Sec. 3. 

*In case of attachment-detachment kinetics, a slope dependent mobility parameter |^ can be used to derive 
the evolution equations. These equations will be reported elsewhere. 
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If the problem under consideration is the decay of a single Fourier mode of amplitude 
uq, the initial conditions satisfy 

«n(o) = h' V (17) 

^ ' ] otherwise. ^ ^ 

An important point to note is that the evolution equations are well-behaved and do not 
directly inherit any singularities that are present in the chemical potential. However, the 
influence of the physical features that give rise to the singularities are incorporated. 

Insights into the nature of the solutions can be gained by expressing (p!6|) in terms of 
dimensionless variables. Introducing the rescaled variables 

h = kh, dn = kttn, X = kx, t = — t (18) 

vr 

we find the Fourier coefficients satisfy 

an{t) = ^dn[pi,i) (19) 

where Pi = /?i//?3 and the functional form of the solution d„ must be determined by solving 
(p!^). We also point out that when the surface shape is not symmetric about the origin, 
a complete Fourier expansion that involves both cosine and sine terms should be adopted; 
the evolution equations for the expansion coefficients can be derived by following the steps 
outlined in this section. 



2.2 Two-dimensional profiles 

The derivation of evolution equations for two-dimensional surface modulations closely follows 
the procedure used in the one-dimensional case. In the case of symmetric profiles, the surface 
shape can be expressed in terms of the Fourier series 

N N 

h{-K,t) = ^ Amnjt) cos{'mkiXi) cos{nk2X2), (20) 

m=l n=l 

where ki^2) = 27r/Ai(2) is the wavenumber of the modulation in the xi(x2)-direction. If the 
amplitude of the modulation is small, we can invoke conservation of mass, expressed in terms 
of the surface shape as 

^+V.j(x,t) = 0, (21) 
to obtain the components of the surface mass flux 

/• v^v^/^- sin(mA;iXi) cos(nA;2X2) ; , , cos(m/i;iXi) sin(n/c2X2) \ 

(jl,j2) = -2^2^ amn{t) 7 ,bmn{t) 7 > (22) 

rr, n \ TTlKl nK2 J 
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where 

O'mn^t) ~\~ ^mnit) — -^mnii)- (23) 

For crystal surfaces that are below the roughening temperature, the surface energy can 
be written as 

= ml + hi)'/' + ^{hl + hlfl\ (24) 

where /3i and /^s have the same meaning as in (^. If we restrict attention to the terrace- 
diffusion limited surface kinetics, the mobility parameter in (^) becomes independent of surface 
slope. In this case, using ( p2[ ) and (24), the variational functional can be written as 



N N N N ^2 

$[dii,di2,...,6ll,...] = - ^ ^{amn + bmn)Hmn + ^ ^ 2k k C 

n=l m=l n=l m=l ^ ^ 

where 



d^ 62 

"'mn _|_ mn 



, (25) 



Hmn 



/ — sm.{mkixi) cos{nk2X2) + nk2-^ — cos(mA;iXi) sin(nA;2a;2) 

V ohx^ dhx2 



dxidx2- 



(26) 

Minimizing the functional with respect to amn and bmm we obtain the first-order ordi- 
nary differential equations 



Adding these two contributions according to (23), we find that the evolution equation for the 
Fourier amplitude in (^) is 

■ C{m'^kfk2 + n'^kik^) 

■^mni^) 9 -f^mn- (28) 

If the amplitudes of the Fourier modes at t = are known, ( |2^ ) can be integrated to obtain 
the surface shape for t > 0. 

The dependence of the decay rate on the aspect ratio of the profile a = X1/X2 can be 
understood by using the following dimensionless variables: 

^1(2) — ^1(2)^1(2); h = kih, Amn = kiA^ani t = 1-, Hmn = -j^Hmn- (29) 

The evolution equations for the scaled Fourier amplitudes in terms of these rescaled variables 
are 

dAyyj.n 



dt 



{m' + a'n')Hmn, (30) 
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where 



Hmn 



2- r2n0^ + hl+^J^ 



mhx^ sin(mxi) cos(nx2) + na^hx^ cos(mxi) sin(nx2) 



dx^ dxo 



(31) 

As noted earlier, the profile decay for the case of /3i = is qualitatively different from the case 
when this parameter is non-zero but arbitrarily small. It is also of interest to examine the 
limiting behavior of the decay of the two-dimensional profiles as a — > 0, which corresponds 
to the situation of one-dimensional behavior determined as a limit of two-dimensional behav- 
ior. To this end, we proceed to show that the two-dimensional evolution equations precisely 
approach the one-dimensional equations as q ^ 0. 

2.3 Decay of elongated 2D profiles: approach of the ID limit 

The approach to the one-dimensional limit can be studied by looking at the decay behavior 
of /i(xi,0,t), given by|] 

h{xuO,t) =J2(fl ^mnityj cos(mxi) (32) 

m \n=l / 

when a ^ 0; the goal is to establish that the quantity A^(i) = J2m=i ^nm(i) obeys the same 
evolution equations as the one-dimensional expansion coefficients, an{t) in (|l^). It follows 
from ( pO| ) and (^) that when a — >■ 0, the scaled Fourier amplitudes satisfy 

"^-m'H^, (33) 



dt 



where 



= m'S^ / / {l3i + h\ )'5gn[hxi]sm{mxi)cos{nx2)dxidx2. (34) 
The integral in ( p^ ) can be explicitly evaluated by noting that 

Sgn[^iJ = Sgn[/i2^]Sgn[cos(i2)], (35) 

where = — X]pP^p° sin(p5:i). In what follows we will separately calculate the contribution 
to due to terms that arise from step formation and interactions. 

^ Since the variational approach is strictly exact only when all the Fourier modes are retained, the upper 
limit in the sum over the Fourier modes will be taken to be infinity. This allows us to use closed form results 
for certain infinite series; these results are required for establishing the limiting ID behavior. 
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Using the results 



27r 



Sgn[cos(x2)] cos{nx2)dx2 



— n = 1, 5, 9, 

n 



n 



n = 3,7, 11, • 
n = 2,4,6, ••■ 



(36) 



and 



E 



(-ir 



vr 



-(2n + l) 4' 

the contribution to (|34|) due to the step formation term is shown to be 



rj-OO I 

\f 



7rm/?i / Sgn[/i|°] sin(mxi)(ixi. 
Jo 



(37) 



(38) 



The contributions due to step interactions can be evaluated by employing the expansion 

^ii ~ E M^pr^qs sin(pxi) sin((7Xi) cos(rx2) cos(sx2) (39) 



p,q,r,s 



and the result 



°0 n2TT 

/ Sgn[cos(x2)] cos(/x2) cos(mx2) cos{nx2)dx2 = tt. 



(40) 



which can be proved for any integers m and n by invoking certain trigonometric identities. 
Using the fact that the sum in (|40| ) is independent of both m and n, it can be shown that the 
contribution to (|3^ due to the step interactions is 



rVco I 



TTm I (/i|^) Sgn[/i2^] sin(mxi)dxi. 



Combining (38) and (|4T|) according to (^4|) , we find that 



TT 



+ Sgn[/i~]sin(nxi)(ixi, 



(41) 



(42) 



which shows that the coefficients A"^ satisfy the same evolution equations as the 1-dimensional 
expansion coefficients (refer to ([T6|)) if we identify Ci, ki and Ai with C, k and A respectively. 
Thus, in the limit of vanishing wavelength aspect ratio, the two-dimensional case reduces to 
the independently established one-dimensional case. Numerical results for profile decay as a 
function of the aspect ratio are discussed in the following section. 
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2.4 Sinusoidal profiles on vicinal surfaces 



As noted earlier, experiments on ID sinusoidal modulations are always affected by the presence 
of steps due to small miscuts in directions in a direction perpendicular to the direction of the 
surface modulations. In this case of a modulated vicinal surface, the surface height can be 
expressed as 

N 

h{xi,X2) = ^ Unit) cos{nkxi) - soX2, (43) 

n=l 

where tan~^(so) is the miscut angle. At time t = 0, we assume that all the Fourier coefficients 
except the longest wavelength mode (n = l)have vanishing amplitudes. A schematic of this 
initial surface shape is given in Fig. |. We win now present the evolution equations for the 
Fourier coefficients and discuss their behavior as sq ^0. 

Using the procedure used to derive the evolution equations for the 2D profiles, we find 
the Fourier coefficients satisfy the evolution equations 

an(t) = ^ ^ Hn[ai, . . . , oat], (44) 

TT 



/■A ^ 

Hn[ai, . . . ,aN] = nk ={j3i + Pshl ) sm{nkxi)dxi. (45) 

^° J hi + si 



where 



It can be seen that these equations have the same form as the evolution equations of the 
Fourier coefficients for the ID profiles given by (ll6|), except for the fact that the Sgn function 
of the surface slope in the expression for Hn is replaced by an analytic function in (|45|). 



Furthermore, since 



hni ^ii== = Sgn[/i,J, (46) 



the evolution equations for the Fourier coefficients in (^) approach the corresponding equa- 
tions for the for the ID modulations derived earlier. We therefore expect the profile shapes on 
the miscut surfaces to coincide with the ID decay profiles of the as the angle of miscut van- 
ishes. The details of the numerical calculations that allow us to look at this limiting behavior 
are given in the following section. 



3 Numerical Results 

In this section, we present results on the decay toward morphological equilibrium of surface 
modulations obtained by numerically integrating the evolution equations derived in Sec. 2. 
We will first consider the decay of one dimensional sinusoidal perturbations for different values 
of step-formation and step-interaction energies. Numerical studies on the sensitivity of the 
results to the number of Fourier components that are used to express the surface shape will 
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then be presented. This is followed by an example of the decay of an anisotropic profile where 
the location of the facets are not easily determined at the beginning of the calculation. In 
the particular case that we consider, some of the facets that appear during the early stages 
of profile evolution, eventually disappear as other neighboring facets grow at faster rates. 
Next, we consider the decay rates of two-dimensional profiles with different aspect ratios and 
diffusion constants. Finally, we focus on the evolution of sinusoidal modulations on a vicinal 
surface that is slightly misoriented with respect to the flat singular surface. 



3.1 ID sinusoidal profiles 

The evolution equations given by (^) were integrated using a Runge-Kutta integrator with 
adaptive step-size control discussed in detail in Ref. [33|. The time dependence of the am- 
plitude decay, for different values of rescaled step-formation energies is shown in Fig. ||. 
The profile shapes at different stages of morphological equilibration are shown, for certain 
representative values of f3i, in Fig. |^. When the step formation energy vanishes, we find 
that the profile sharpens at its extrema and the decay rate shows the 1/t scaling behavior 
predicted by Ozdemir and Zangwill ||l^ (refer to inset (a) in Fig. ^). On the other hand, for 
any non-zero value of we find that the profile decays with the formation of facets which, 
after an initial transient period of growth, remain nearly constant in size as seen in Fig. 
The amplitude during this self-similar decay stage, decreases as a linear function of time, 
(refer to Fig. |2|), in agreement with the results of Hager and Spohn which were obtained 
using the free-boundary approach. With increasing value of (5i, the lengths of the facets in 
the self-similar shape increase in magnitude and the initial transient periods decrease in their 
duration. These results indicate that the case /9i 7^ is a singular perturbation of the special 
case with (3i = 0; any non-zero value of /3i leads to facet formation during the morphological 
equilibration. 

The sensitivity of the profile shapes and decay rates on the number of Fourier modes 
used in the expansion of the surface shape (refer to (|l2[)) is considered in Fig. ^, where the 
decay rates obtained by retaining 15, 20 and 30 modes are given. The insets in the figure 
show the evolution of the amplitudes of the individual Fourier modes for the three cases 
considered. We find that the amplitude decay remains virtually unchanged as the number of 
Fourier modes increase from 15 to 30. Fig. |^ shows that the dominant Fourier components 
in the surface shape are the long wavelength modes and that the magnitudes of the different 
Fourier components are also seen to be independent of the total number of modes that are 
used in the analysis. The insets also show that all the Fourier components decay linearly in 
time after an initial transient period; in this regime, the surface shape remains self-similar 
as noted earlier. We have investigated the time dependence of the decay shapes for several 
other values of the step-formation energies and initial amplitudes of the sinusoidal profiles 
and found that accurate results are obtained by using between 15 and 20 Fourier modes. As 
is the case with any type of spectral method, if more basis functions are retained, the time 
steps used in numerical integration have to be decreased, which leads to longer simulation 
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times. The optimum choice for a particular application is determined by a trade-off between 
accuracy and simulation time. 

3.2 Asymmetric profiles 

Next, we consider the evolution of an asymmetric profile shown in Fig. ^(a), where the surface 
shape can be expressed by using both cosine and sine functions in the Fourier expansion as 

h{x, t) = ^ am{t) cos{nkx) + ^ bn{t) sm{nkx). (47) 

m n 

The amplitudes of the non-vanishing Fourier components in the initial shape were chosen as 
follows: kai = 0.2, ka^ = 0.05, kb2 = 0.1 and kb^ = 0.025. We find that during the early stages 
of morphological evolution, facets labeled 1, 2 and 3 appear as shown in Fig. ||(b). While these 
facets are situated close to the extrema in the initial profile, their actual locations are not 
determined very easily. As the shape evolves in time, the facet labeled 3 shrinks in size, while 
the other two facets grow is size; in Fig. ^(c) this facet has almost disappeared. In the late 
stages of the evolution, all the Fourier sine components with short wavelengths vanish, so that 
the shape closely resembles the decay profile of the symmetric cosine mode with the largest 
wavelength. 

The above example illustrates that the present approach to surface evolution can be 
used to study arbitrary profile shapes in a straightforward manner. While the free-boundary 
approach, in principle, can be used in this case, the actual implementation can be quite tedious. 
First, the locations of the facets are difficult to guess at the outset and secondly, application of 
boundary conditions at the points at which the facets disappear require special consideration. 
On the other hand, there is no fundamental difference in the analysis of symmetric and 
asymmetric profiles within the variational approach. 

3.3 2D sinusoidal profiles 

In the case of bidirectional modulations, we start by looking at the decay behavior of a 
bidirectional profile whose wavelengths of modulation in the coordinate directions are equal. 
The amplitude evolution and profile shapes for different values of f3i are given in Fig. ^ and 
Fig. ^, respectively. Just as in the ID case, when /3i = 0, we find that the profile decays 
without formation of facets. Furthermore, the amplitude of the modulation decays with the 
1/t scaling behavior, which is also observed in the decay of ID profiles. It is also evident from 
Fig. that the surface profile close to extrema becomes shaper than the parabolic shape, 
which once again resembles the behavior in ID. For any non-zero the profile decays with 
the formation of facets, with the decay rate increasing with Pi. The decay rate, however, does 
not follow the 1/t scaling behavior; instead, it shows a linear behavior after an initial transient 
as indicated in Fig. 0; this trend is also similar to the decay of a ID profile. These results 
indicate that there is no fundamental difference in the profile decay behavior in one- and two- 
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dimensions; the decay behavior is governed primarily by the ratio of the step formation and 
interaction energies. 

In Sec. 2, we presented analytical arguments that establish that the decay of elongated 
the 2D profiles approaches the ID limit as the aspect ratio of the profile diverges. We now 
discuss numerical calculations that allow us to understand the way the profile shapes evolve 
with increasing aspect ratios. In Fig. the decay rates for profiles with different values 
of a are given along with the curve for the case a = 0, calculated with the ID evolution 
equations (refer to (16)). We find that as a falls below 1, the decay rates initially decrease 



continuously until a reaches about 0.001. As a decreases further, there is no significant change 
in the profile decay, so that in the limit of very small a all the decay curves converge to a 
single curve. It can be seem from Fig. ^ that this curve coincides with the decay profile of 
the ID case, in agreement with the analytical results of Sec. 2. We have also looked at the 
convergence properties for several other values of f3i and initial amplitudes and have found 
that the ID limit is achieved when a lies between 0.01 and 0.001. The evolution of the surface 
shape during the decay of an elongated profile with a = 0.1 is illustrated in Fig. |8|. It can 
be seen that well-defined straight features that are aligned along the direction in with the 
profile is elongated, appear in the decay shape. The cross-section of the profile, particularly 
at the extrema of the elongated shape, are found to closely resemble the decay profile of a ID 
modulation. 

3.4 Sinusoidal profiles on vicinal surfaces 

As a final example of the application of the variational approach, we consider the decay 
of sinusoidal profiles on vicinal surfaces that are slightly misoriented from the fiat singular 
surface. Our goal is to understand the limiting behavior of profile decay as the miscut slope 
of the surface approaches zero. To this end, we consider the decay of a sinusoidal modulation 
with scaled amplitude, kai = 0.1 in Fig. The step formation energy was taken to be 
/?! = 1.0. The figure shows decay curves for several values of the misorientation slopes, along 
with the decay profile for the ID modulation calculated using (p!6|). We find that the profile 
decay progresses at a faster rate as the misorientation slope decreases from 0.1 to 0.01. As 
the slope is decreased further, all the decay curves approach a limiting curve which coincides 
with the ID decay profile. In the particular example that we consider, the decay behavior 
is already very close to the ID case when the misorientation slope is about 0.001. We have 
also numerically investigated the decay behavior for several other values of step formation 
energies and initial amplitudes and found that the decay curves approach the ID case as the 
misorientation slope approaches zero. 

Our results differ from the earlier results of Lancon and Villain [|l^], who point out 
that facets form on surfaces only if the misorientation slope is larger than some critical value. 



For the range of parameters considered in Fig. |10|, they estimate the critical value to be 
■So = 0.001. In contrast to their predictions, our numerical studies clearly show that facets 
form for arbitrarily small values of the miscut slope. This observation is consistent with 
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the fact that the evolution equations for the Fourier coefficients that determine the shape 
modulations on miscut surfaces approach the corresponding equations for ID coefficients as 
the miscut slope vanishes. 



3.5 Comparison to experiments on Si(OOl) 

We now turn our attention to recent experimental work of Erlebacher and coworkers |31|, who 
investigated the decay of sputter ripples on Si(OOl) surfaces. They found that amplitude of 
the ripples, with a = A1/A2 = 0.1, exhibited a decay behavior consistent with the 1/t scaling 
behavior for ID modulations predicted by Zangwill and Villain and cowor kers ||,||. In this 
section, using available experimental and theoretical data on step formation and interaction 
energies for Si(OOl), we will numerically investigate the decay of these elongated ripples. 

The parameters /3i and /?3, for the Si(OOl) surfaces with alternating SA and SB steps 
relevant to the experiments, can be estimated from experimentally determined step formation 
energies reported by Zandvliet [34| and the atomistic calculations of Poon and coworkers [Q. 
Using the step formation energies of SA (8 me V/A) and SB (16 meV/A) steps, we find that 
/?! sa 12.5 meV/A^. From the atomistic calculations of step interactions on Si(OOl) |3^, /^s is 
estimated to be ~1200 meV/A^. The profile decay for the elongated ripples with /3i = 0.01 is 



shown in Fig. 11. As noted earlier, when /3i = 0, both ID and 2D profiles show a 1/t decay 
behavior. The figure shows that in distinct contrast to the /3i = case, even for a relatively 
small non-zero value of the amplitude decay does not exhibit the 1/t scaling behavior. The 
observed behavior is therefore not consistent with diffusion limited decay of 2D ripples. The 
results can perhaps be explained by including anisotropics in surface diffusion and the effects 
of force-monopoles at the steps due to the 2x1 and 1x2 reconstruction of the terraces between 
the steps. It has also been suggested that the surface transport at the temperatures of interest 



is limited by attachment-detachment kinetics [36|. We leave a detailed consideration of these 
issues for future study. 



4 Concluding Observations 

In conclusion, we have developed a numerical method to study the singular equations that 
govern the dynamics of crystal surfaces below the roughening temperature. The method 
was applied to study the morphological equilibration of both unidirectional and bidirectional 
modulations. In all these cases, the surface evolution was determined by solving well-behaved 
coupled nonlinear ordinary differential equations for the modal expansion coefficients of the 
surface shape. Analysis of highly elongated 2D sinusoidal profiles and sinusoidal profiles on 
vicinal surfaces with vanishing miscut angles shows that the limiting behavior in both of these 
cases agrees with the ID decay behavior accompanied by facet formation. In the following 
paragraph, we point to directions of further work that will allow us to make closer connection 
to experiments. 
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In this paper we restricted attention to surface diffusion kinetics limited by terrace 
diffusion, while many experimental systems are in a regime where the mass transport is 
limited by attachment and detachment at the steps [|lO|. For the latter case, Nozieres has 
derived a kinetic relation between the gradient of the surface chemical potential and the mass 
flux . In the future we plan to derive a variational formulation appropriate to this kinetic 
law to analyze morphological evolution in the attachment-detachment limited regime. Also, in 
some experimental situations, anisotropics associated with surface diffusion as well as surface 
energies become can be significant. As noted earlier, these effects can be included in the 
present formulation and will be the subject of forthcoming publications. 
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FIGURES 




Figure 1: Sketch of a sinusoidal profile on a surface that is slightly misoriented from 
the flat singular surface. The wavelength of the perturbation is A and the spacing 
between the steps that arise due to miscut is d = hg/so, where hs is the step height 
and So is the slope of the miscut surface relative to the flat orientation. 
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Figure 2: Normalized amplitude h{0,t)/ao as a function of rescaled time t for the 
decay of ID sinusoidal profiles plotted for different values of step-formation energies, 
Pi- The initial amplitude of the profile was taken to be a(0) = 0.1 /k, where k is 
the wavenumber. The inset labeled (a) shows the time dependence of both the 
normalized amplitude and its inverse for the case of zero step-formation energy 
(/?! = 0). Here, amplitude decay is seen to scale as h{0,t) ~ 1/t, in agreement with 
the results of |12| and |14|. Inset (b) shows that this behavior does not, however, 
hold for non-zero values of f3i. The dotted lines next to the decay curves indicate 
that the amplitudes, after an initial transient period, decay linearly in time. 
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Figure 3: Morphological equilibration of a symmetric sinusoidal profile for different 
values of the step- formation energy, Note that the propensity to form flat tops 
increases with increasing values of Pi. When Pi =0.1 and 1.0, the shapes of the 
curves 1 and 2 marked in the figure are self-similar; the lengths of the facets remains 
constant during this decay stage. 
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Figure 4: Decay of a sinusoidal profile calculated by using 15, 20 and 30 Fourier 
modes to express the surface shape. The initial amplitude and the step formation 
energy(in rescaled units) were taken to be ai(0) = 0.1 and /3i = 1.0, respectively. 
The figure shows that the profile decay is virtually unchanged as the number of 
Fourier modes are doubled, indicating convergence of the method can be achieved by 
retaining 15-20 modes. The insets in the figure show the evolution of the amplitudes 
first 6 non-zero odd modes for the three cases considered (symmetry arguments can 
be invoked to show that even modes do not grow in amplitude). The mode labels 
are given in the inset corresponding to N=15. Note that the dominant modes in the 
profile shapes are the ones whose wavelengths are large. Also, the amplitudes of the 
individual modes remain unchanged as the total number of modes increases from 15 
to 30. 
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Figure 5: The sequence of frames show the time evolution of an asymmetric profile 
that contains both cosine and sine modes. The initial profile shape is shown in (a); 
the amplitudes of the non-vanishing modes at t = are kai = 0.2, ka^ = 0.05, kb2 = 
0.1 and kb^ = 0.025 (refer to (p7[)), where k is the wavenumber of the perturbation. 
During early stages of decay, three facets labeled 1, 2 and 3 appear as shown in (b). 
As the evolution proceeds, the facets labeled 1 and 2 grow at the expense of 3, which 
has almost disappeared in (c) . In the late stages of the evolution shown in (d) , the 
shape closely resembles the decay profile of the symmetric longest wavelength cosine 
mode. In this calculation, we have chosen /3i = 5.0. The profiles in (a),(b),(c) and 
(d) are plotted at i = 0.0, 1.9 x 10-^ 9.5 x lO'^^ and 3 x 10-^ respectively. 
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Figure 6: Normalized amplitude h{0,t)/aQ as a function of rescaled time t for the 
decay of 2D sinusoidal profiles plotted for different values of step-formation energies, 
$1 . The wavelengths of the profiles in the coordinate directions are taken to be equal 
and initial amplitude was assumed to be ao = O.l/fc, where k is the wavenumber. 
The inset labeled (a) shows the time dependence of both the normalized amplitude 
and its inverse for the case of zero step-formation energy = 0). Here, just as in 
the ID case, the amplitude decay is seen to scale as h{0,t) ~ Inset (b) shows 
that this behavior does not, however, hold for non-zero values of The dotted 
lines next to the decay curves indicate that the amplitudes, after an initial transient 
period, decay linearly in time. 
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Figure 7: Time sequence of the decay of bidirectional sinusoidal profiles with a = 
A1/A2 = 1, where A, are the modulating wavelengths in the coordinate directions. 
The initial amplitude of the modulation is chosen such that kAn = 0.1, where k is 
the wavenumbcr of the profile. When (3i = 0, the profile decays without formation 
of any facets, while facets are formed for any non-zero value of For the case 
$1 = 1.0, formation of facets at the extrema of the profile is clearly seen in the 
figure. 
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Figure 8: Time sequence of the decay of a two-dimensional profile with A1/A2 = 0.1, 
where Aj are the modulating wavelengths in the coordinate directions. The scaled 
step- formation energy is taken to be Pi = 1.0 and the initial amplitude is chosen 
such that kiAii = 0.1, where ki is the wavenumbcr along the 1-dircction. As the 
profile evolves, elongated ridges and troughs, aligned along the 2-direction appear in 
the surface shape. The cross-section of the profile, particularly at the extrema of the 
elongated shape, is found to closely resemble the decay profile of a ID modulation. 
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Figure 9: Profile decay for bidirectional sinusoidal modulations plotted for different 
values of the profile aspect ratios. The aspect ratio is varied by holding the wave- 
length in the xi-direction fixed while increasing the wavelength of modulation in 
the X2-direction. The scaled step-formation energy is taken to be Pi = 1.0 and the 
initial amplitude is chosen such that kiAn =0.1, where ki is the wavenumber along 
the 1-direction. The figure also shows profile decay for a ID sinusoidal modulation 



(dashed-line) calculated using (16). It can be seen that the decay curve of the highly 
elongated profile with a = 0.001 very closely follows that for the ID case. 
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Figure 10: Profile decay for sinusoidal modulations on surfaces that are slightly 
miscut from the flat orientation, plotted for different values of miscut slopes (sq). 
The initial amplitude and the step formation energy (in rescaled units) were taken 
to be ai(0) = 0.1 and f3i = 1.0, respectively. The figure also shows profile decay 
for a ID sinusoidal modulation (sq = 0) calculated using (1€). Note that the decay 
curve of the vicinal surface with sq = 0.001 is seen to almost coincide with the ID 
curve. The decay curves for sq < 0.0001 have not been displayed since they almost 
exactly coincide with the ID decay curve. The insets in the figure show the level 
curves of surface height (refer to (|43|)) when sq = 0.001. The inset labeled (a) shows 
the level curves of the starting sinusoidal profile, while (b) shows the level curves 
after the formation of facets at the extrema of the profile. 
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Figure 11: Profile decay of sinusoidal modulations with A2/A1 = 10, considered 
in the experimental work of Erlebacher and coworkers ||3l|. The figure also shows 
the time dependence of the inverse amplitude for two different values of the scaled 
step formation energies. In both the cases, the amplitude does not show at the 1/t 
behavior observed in the experiments. Using experimentally determined vales of 
step formation energies [34| and step interaction energies obtained form atomistic 
calculations we find that /3i « 0.01 for SA+SB steps on Si(OOl). 
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